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Abstract 

Recent experiments by Baxter et al. showed the existence of density waves in granular material 
flowing out of a hopper. We show, using Molecular Dynamics Simulations, that this effect is a 
consequence of static friction and find that these density fluctuations follow a 1/ f spectrum. The 
effect is enhanced when the opening angle of the hopper decreases. 
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Moving dry granular media, like sand, show a rich variety of rather astonishing and scarcely 
understood phenomena' 1-2 !. Famous are the so-called "Brazil nut" segregation! 3-5 ! and the heap 
formations that occur under vibrations' 6-8 !. More recently a series of experiments have given evi- 
dence that under certain circumstances density patterns are generated inside the flowing medium. 
Baxter et alJ 9 ', for instance, visualized wave-like patterns emanating from the outlet of a two 
dimensional wedge-shaped hopper using X-rays. Also previous authors' 10-12 ! had noted the for- 
mation of similar structures. Similarly rather erratic shock-like density waves have been observed 
in flow through pipes' 13 ! and down inclined planes' 14 ! . Another experimentally observed ubiquous 
phenomenon in granular media seems to be l/,f a noise. Lui and Jaeger' 15 ! recently measured the 
acceleration of a particle inside a bulk of glass beads that where excited by a small amplitude 
vibration. Its Fourier spectrum in time showed power law decay over many orders of magnitude. 
Baxter' 16 ! also observed power law decay in the frequency dependent forces that act on the wall 
of a hopper. For avalanches going down the slope of a sand pile theoretical considerations of self- 
organized criticality' 17 ! led to the proposal that their size and life times distributions were power 
laws which was in fact only verified experimentally on very small piles' 18 ' . 

The existence of these erratic density and force inhomogencities arc intimately related to the 
ability of granular materials to form a hybrid state between a fluid and a solid: When the density 
exceeds a certain value, the critical dilatancy' 19 ' 20 ! , granular materials are resistant to shear, like 
solids, while below this density they will "fluidify". In the presence of density fluctuations the 
rheology therefore can become rather complex. Two microscopic facts seem to be responsible for 
the strong density fluctuations: On one hand one has in granular media solid friction between 
the grains. This means that when particles are pushed against each other a finite force is needed 
to start or maintain a relative tangential motion between them. On the other hand a granular 
material is internally disordered giving a natural source of noise. The strong non-linearity coming 
from friction produces instabilities in the density which enhance the fluctuations coming from the 
noise. 

Various attempts have been made to formalize and quantify the complicated rheology of gran- 
ular media. Continuum equations of motion ' 21 1, a cellular automaton ' 22 1 and a random walk 
approach' 23 ! have been proposed. But none of them has yet been able to explain these heteroge- 
neous waves. This is why we chose to study these phenomena using Molecular Dynamics (MD) 
simulations of inelastic particles with static and dynamic friction in two dimensional systems. In 
fact, MD simulations' 24 ' 25 ! have already been applied to granular media to model segregation' 5 !, 
outflow from a hopper' 26, 27 1, shear flow' 28 !, convection cells on vibrating plates' 29 ' 30 !, avalanches 
on a sand pile' 31 ! j fl ow through a pipe' 32 ! and others. 

We consider a system of N spherical particles of equal density and with diameters d cither 
all equal or chosen randomly from a Gaussian distribution of width w around do = 1 mm. These 
particles are placed into a hopper having an opening angle 9 and at the bottom an opening of 
diameter D. When two particles i and j overlap (i.e. when their distance is smaller than the sum 
of their radii) three forces act on particle i: 1.) an elastic restoration force 

f1?=Y(\n j \-\(d i + d j )) J ^- , (la) 

z I >ij I 
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where Y is the Young modulus and points from particle i to j; 2.) a dissipation due to the 
inelasticity of the collision 

fdL = -7 m eff(Vij ■ T-T-H = -7 m effVij , (16) 

I 'ij I 

where 7 is a phenomenological dissipation coefficient and wjj — Vi~ Vj the relative velocity between 
the particles; 3.) a shear friction force which in its simplest from can be chosen as 

flhear = -1*™>eff tfij ' Uj ) l^TTi = -1> m eff V ij ' ( 2fl ) 

I 'ij I 

where 7 S is the shear friction coefficient and tij = (— r\j,rfj) is the vector r*jj rotated by 90°. 
Eq. 2a is a rather simplistic description of shear friction. In our case it is important to include 
real static friction which is done by a static friction force' 33 ': When two particles start to touch 
each other, one puts a "virtual" spring between the contact points of the two particles. Be 6s the 
total shear displacement of this spring during the contact and k s 5s the restoring frictional force 
(static friction). The maximum value of the restoring force is then according to Coulomb's criterion 
proportional to the normal force F n and the proportionality constant is the friction coefficient /i. 
Cast into a formula this gives a friction force 

ffriction = -sign(6s)mm(k s m eff 6s, ^F n ) . (26) 

where 6s is the shear displacement integrated over the entire collision time. When particles are 
no longer in contact with each other the spring is removed. It is, however, not straightforward to 
implement the above technique when the particles are allowed to rotate, i.e. to be able to roll on 
each other. Therefore we did not take angular momenta into account. In fact, when particles have 
strong deviations from the spherical shape rotations are strongly suppressed. Often it is however 
useful to go an intermediate way and to include dynamic friction but not static friction and allow 
for the particles to have rotations' 5,27,29 '. In that case one uses a combination of eqs. 2a and 2b: 

fdyn = -min(7 s m e// 4 , fiF n ) (2c) 

and introduces also equations of motion for the angular momentum of the particles. 

When a particle collides with a wall the same forces act as if it would have encountered another 
particle of diameter do with infinite mass at the collision point. The walls are in fact made out 
of small particles themselves and in order to introduce roughness on the wall these particles are 
chosen randomly from a distribution of two radii. The only external force acting on the system is 
gravity g ~ — 10m/s 2 . 

As initial positions of the particles we considered that they are placed at random positions 
inside a space several times as high as the dense packing. The initial velocities are set to zero. 
After that the particles are allowed to fall freely under gravity. We took a Young modulus of 
Y = 10 6 g/s 2 and a time step ofAi = 3.3-10~ 6 s and run our program on 8 or 16 processors of an 
Intel iPSC/860 and an IBM RS/6000-550. 

In Fig. 1 we see a snapshot of the outflowing particles at four different time steps for material 
parameters consistent with the experiment of Baxter et al.' 2 ' . We clearly see that close to the outlet 
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large holes appear which then propagate in an attenuated form upwards. These patterns quickly 
vary in time. When the opening angle 9 of the hopper is reduced the contrast in the patterns 
becomes more pronounced. The holes are in general stretched in the horizontal direction but their 
shape seems rather random. Particularly striking are these structures when watched in a movie. 
When the static friction \i is switched off the structures disappear and the density of the outflowing 
particles becomes homogeneous. This agrees with Baxter et al.'s experimental observation that 
density patterns only occur for rough and not for smooth sand. When the friction with walls is 
switched off the density waves also disappear. 

A more quantitive approach can be made by measuring the local densities p. We binned space 
in units of 1.56do and counted the number of particles that have their center of mass inside the 
box averaged over 100 consecutive iteration steps. This density is plotted in Fig. 2 as a function of 
time and space. We see that the low density regions form curved stripes pointing to rather short 
lived waves. No regular structure in the distance and magnitude of the waves can be noted. They 
rather look like independent shock waves with random amplitudes coming in a random sequence. 
Larger waves usually have small densely packed precursors. Strong similarities can in fact be seen 
with space-time plots of the density of granular media in pipes I 32 l and of cars on highways having 
traffic jams! 34 ! . 

In Fig. 3a we see the density profile as a function of time for a point four particle diameters 
above the outlet. No regularity can be seen. The over four components averaged Fourier trans- 
formed data are shown in a log-log plot in Fig. 3b. Clearly they fall on a straight line over nearly 
two decades. The slope is about -1. This means that we have found 1// noise. This is consistent 
with the picture that has emerged from the mentioned experiments. 

When particles of equal size are taken we observed equally well developed density patterns and 
find roughly the same power law decay of the spectrum. The effect is reduced when the diameter 
D of the outlet becomes too large. If it is too small the flow of sand can entirely stop due to 
arching. The critical diameter Dq when this arching sets in has been studied before with similar 
techniques I 26 ' where it was found that Do is larger when the particles have the same size. When 
we consider smooth walls, i.e. all wall particles having the same radii, we do not find density 
waves and the power spectrum looks significantly different. It shows an upwards curved slope with 
increasing frequency which one also finds when configurations block during the outflow. A similar 
effect was also found in simulations of flow on an inclined plane! 35 ! _ 

Baxter ct al. also observed that there exists a critical angle 9q above which stagnation regions 
exist next to the walls of the container. We also found these regions one example is shown in Fig. 4. 

We have observed in a very simple modclization that density waves are generated and dis- 
tributed with frequency like 1//. Two ingredients were found essential to generate them: static 
friction and disorder. The static friction tends to align the particles, i.e. to form fronts of particles 
moving exactly with the same vertical velocity. These fronts are nucleated randomly at the walls. 
Their size distribution (density contrast) comes by itself into a critical state namely a power law 
distribution. It therefore has the properties of self-organized criticality (SOC)! 17 '. It is, however, 
very important to notice that our simulations were made for rather small systems as compared to 
real systems. It could therefore be that for systems of millions of particles a cut off exists in this 
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power law. 

We have shown in this paper that similar to the avalanches that one observes on the surface 
of a sand pile also inside the bulk of granular material one has avalanche behavior which like the 
ones on the surface shows self-organized criticality on small scales f 18 '. Surface avalanches for large 
piles however seemed to have a characteristic size due to inertia! 36 !. j n f ac t ; rather avalanches 
in the holes than in the mattered substance could be relevant here and one might speculate that 
the bulk avalanches might be a better example for asymptotic SOC than the ones on the surface. 
The mechanisms that generate the patterns are similar but not identical to the original sandpile 
models. While the static friction similarly generates waiting times with a threshold it is not the 
motion of the sand itself that constitutes the avalanches but it is the group velocity of the holes 
between them: An individual particle can easily go from one dense region to the other by flying 
fast through a region of low density. There is therefore a backflow of information similar to the 
jamming on highways I 34 '. 

Although our simulations are two dimensional, we think that they do capture the essential 
mechanisms that occur also in three dimensional experiments. It should be mentioned that in 
fluidized beds (low Bagnolds number) where the granular medium is surrounded by a fluid and 
the hydrodynamic interactions become important a similar phenomenon as the one described here, 
called slugging, is observed! 37 ! . The mechanisms involved seem, however, quite different. 



* present address: Groupe de Matiere Condenscc et Materiaux, Universite dc Rcnncs I, Cam- 
pus dc Beaulieu, 35000 Rennes, France 
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Figure Captions 

Figure 1. In this sequence, we show the outflow behavior with static friction for an initial configuration 
of roughly 1500 particles. The parameters were chosen as following: 7 = 500Hz, 7 S = 250Hz, 
k s = 1000g/s 2 ,6> = 30°, D = 8do,M = 0-5. (a) initial configuration, (b) after 10,000 iterations, 
(c) after 100,000 iterations, (d) after 140,000 iterations 

Figure 2. This figure shows the spatial density fluctuations (vertical axis) for a typical outflow simulation 
as a function of time (horizontal axis) for 7 = 100Hz, "f s = 500Hz, k s = 1000g/s 2 , 9 = 10°, D = 
10do and \i = 0.5. Dark areas mark regions with lower densities. 

Figure 3. We consider an outflow sequence for the same parameter values as in figure 1 but we constantly 
fill in particles from above, a) shows the density fluctuations 6c?o above the hole. One clearly 
sees a quite random behavior, b) shows the log-log plot of the FFT analysis of the density 
fluctuations. A straight line with slope -1 is drawn to guide the eye. 

Figure 4. Due to static friction, 314 particles remain in the hopper for an opening angle of 9 = 150° 
(7 = 100Hz, 7s = 0Hz, k s = 1000g/s 2 , D = 10d and = 0.5). 
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